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Abstract 

Point pattern matching problems are of fundamental importance in various areas including computer 
vision and structural bioinformatics. In this paper, we study one of the more general problems, known 
as LCP (largest common point set problem): Let P and Q be two point sets in R 3 , and let e > be a 
tolerance parameter, the problem is to find a rigid motion fj, that maximizes the cardinality of subset 7 
of Q, such that the Hausdorff distance dist(P, /z(7)) < e. We denote the size of the optimal solution to 
the above problem by LCP(P, Q). The problem is called exact-LCP for e = 0, and tolerant-LCP when 
e > and the minimum interpoint distance is greater than 2e. A /3-distance-approximation algorithm for 
tolerant-LCP finds a subset ICQ such that |7| > LCP(P, Q) and dist(P, fi(I)) < (3e for some /3 > 1. 

This paper has three main contributions. (1) We introduce a new algorithm, called DlHEDA , which 
gives the fastest known deterministic 4-distance-approximation algorithm for tolerant-LCP. (2) For the 
exact-LCP, when the matched set is required to be large, we give a simple sampling strategy that improves 
the running times of all known deterministic algorithms, yielding the fastest known deterministic algo- 
rithm for this problem. (3) We use expander graphs to speed-up the DlHEDA algorithm for tolerant-LCP 
when the size of the matched set is required to be large, at the expense of approximation in the matched 
set size. Our algorithms also work when the transformation /i is allowed to be scaling transformation. 

Keywords. Point Pattern Matching, Largest Common Point Set 

1 Introduction 

The general problem of finding large similar common substructures in two point sets arises in many areas 
ranging from computer vision to structural bioinformatics. In this paper, we study one of the more general 
problems, known as the largest common point set problem (LCP), which has several variants to be discussed 
below. 



Problem Statement. Given two point sets in M 3 , P = {pi, . . . ,p m } and Q = {qi, . . . , q n }, and an error 
parameter e > 0, we want to find a rigid motion fi that maximizes the cardinality of subset I Q Q, such 
that dist(P, /x(-T)) < e. For an optimal set 7, denote |7| by LCP(P, Q). There are two commonly used 
distance measures between point sets: Hausdorff distance and bottleneck distance. The Hausdorff distance 
dist(P, Q) between two point sets P and Q is given by max 9g g min^p \ \pq\\- The bottleneck distance 
dist(7', Q) between two point sets P and Q is given by minj max 5S Q \ \f(q) — q\ |, where / : Q — > P is an 
injection. Thus we get two versions of the LCP depending on which distance is used. 

*A preliminary version was presented at the 7th International Symposium, Latin American Theoretical Informatics (LATIN 
2006) [15]. 

' Corresponding author. Department of Computer Science, Virginia Tech, USA. vchoi@cs.vt.edu. 
^Department of Computer Science, McGill University, Canada. navin@cs.mcgill.ca. 
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Another distinction that is made is between the exact-LCP and the threshold-LCP. In the former we have 
e = and in the latter we have e > 0. The exact-LCP is computationally easier than the threshold-LCP; 
however, it is not useful when the data suffers from round-off and sampling errors, and when we wish to 
measure the resemblance between two point sets and do not expect exact matches. These problems are 
better modeled by the threshold-LCP, which turns out to be harder, and various kinds of approximation 
algorithms have been considered for it in the literature (see below). A special kind of threshold-LCP in 
which one assumes that the minimum interpoint distance is greater than the error parameter 2e is called 
tolerant-LCP. tolerant-LCP more accurately captures many problems arising in practice, and it appears that 
it is algorithmically easier than threshold-LCP. Notice that for the tolerant-LCP, the Hausdorff and bottleneck 
distances are essentially the same in the sense that the problem has a solution of Hausdorff distance < e if 
and only if the solution is of bottleneck distance < e. Thus, for the tolerant-LCP, there is no need to specify 
which distance is in use. 

In practice, it is often the case that the size of the solution set / to the LCP is required to be at least 
a certain fraction of the minimum of the sizes of the two point sets: |/| > - min(|P|, \ Q\), where a is 
a positive constant. This version of the LCP is known as the a-LCP. A special case of the LCP which 
requires matching the entire set Q is called Pattern Matching (PM) problem. Again, we have exact-PM, 
threshold-PM, and tolerant-PM versions. 

In this paper, we focus on approximation algorithms for tolerant-LCP and tolerant-ct-LCP. There are 
two natural notions of approximation. (1) Distance approximation: The algorithm finds a transformation 
that brings a set / C Q of size at least LCP(P, Q) within distance e' for some constant e' > e. (2) Size- 
approximation: The algorithm guarantees that |/| > (1 — <5)LCP(P, Q), for constant 5 G [0, 1). 

Previous work. The LCP has been extensively investigated in computer vision (e.g. [31]), computational 
geometry (e.g. [8]), and also finds applications in computational structural biology (e.g. [33]). For the 
exact-LCP problem, there are four simple and popular algorithms: alignment (e.g. [26, 5]), pose clustering 
(e.g. [31]), geometric hashing (e.g. [30]) and generalized Hough transform (GHT) (e.g. [22]). These 
algorithms are often confused with one another in the literature. For convenience of the reader, we include 
brief descriptions of these algorithms in the appendix. Among these four algorithms, the most efficient 
algorithm is GHT. 

Exact algorithms for tolerant-LCP. As we mentioned above, the tolerant-LCP (or more generally, 
threshold-LCP) is a better model of many situations that arise in practice. However, it turns out that it 
is considerably more difficult to solve the tolerant-LCP than the exact-LCP. Intuitively, a fundamental dif- 
ference between the two problems lies in the fact that for the exact-LCP the set of rigid motions, that may 
potentially correspond to the solution, is discrete and can be easily enumerated. Indeed, the algorithms for 
the exact-LCP are all based on the (explicit or implicit) enumeration of rigid motions that can be obtained 
by matching triplets to triplets. On the other hand, for the tolerant-LCP this set is continuous, and hence the 
direct enumeration strategies do not work. Nevertheless, the optimal rigid motions can be characterized by a 
set of high degree polynomial equations as in [9]. A similar characterization was made by Alt and Guibas in 
[7] for the 2D tolerant-PM problem and by the authors in [14] for the 3D tolerant-PM. All known algorithms 
for the threshold-LCP use these characterizations and involve solving systems of high degree equations 
which leads to "numerical instability problem" [7]. Note that exact-LCP and the exact solution for tolerant- 
LCP are two distinct problems. (Readers are cautioned not to confuse these two problems as in Gavrilov et 
al. [18].) Ambiihl et al. [9] gave an algorithm for tolerant-LCP with running time 0(m 16 n 16 ^m + n). The 
algorithm in [14] for threshold-PM can be adapted to solve the tolerant-LCP in 0(ra 6 n 6 (m + n) 2 5 ) time. 
Both algorithms are for bottleneck distances. These algorithms can be modified to solve threshold-LCP 
under Hausdorff distance with a better running time by replacing the maximum bipartite graph matching 
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algorithm which runs in 0(n 2 ' 5 ) with the 0(n log n) time algorithm for nearest neighbor search. Both of 
these algorithms are for the general threshold-LCP, but to the best of our knowledge, these algorithms are 
the only known exact algorithms for the tolerant-LCP also. 

Approximation algorithms for tolerant-LCP. Like threshold-LCP, the exact algorithm for threshold- 
PM is difficult, even in 2D (see [7]). Two types of approximation algorithms were studied. First, 
Goodrich et al [19] showed that there is a small discrete set of rigid motions which contains a rigid mo- 
tion approximating (in distance) the optimal rigid motion for the threshold-PM problem, and thus the 
threshold-PM problem can be solved approximately by an enumeration strategy. Based on this idea and 
the alignment approach of enumerating all possible such discrete rigid motions, Akutsu [4], and Biswas and 
Chakraborty [11, 10] gave distance-approximation algorithms with running time 0(m 4 n 4 ^m + n) for the 
threshold-LCP under bottleneck distance, which can be modified to give 0(m 3 n 4 log m) time algorithm 
for the tolerant-LCP. Second, Heffernan and Schirra [23] introduced approximate decision algorithms to 
approximate the minimum Hausdorff distance between two point sets. Given e > 0, their algorithm an- 
swers correctly (YES/NO) if e is not too close to the optimal value e* (which is the minimum Hausdorff 
distance between the two point sets) and DON'T KNOW if the answer is too close to the optimal value. 
Notice that this approximation framework can not be "similarly" adopted to the LCP problem because in 
the LCP case there are two parameters - size and distance - to be optimized. This appears to be mis- 
taken by Indyk et al. in [25, 18] where their approximation algorithm for tolerant-LCP is not well defined. 
Cardoze and Schulman [12] gave an approximation algorithm (with possible false positives) but the trans- 
formations are restricted to translations for the LCP problem. Given a, let e m j n (a) denote the smallest e 
for which q-LCP exists; given e, let a m i n (e) denote the smallest a for which ct-LCP exists. Biswas and 
Chakraborty [11, 10] combined the idea from Heffernan and Schirra and the algorithm of Akutsu [4] to give 
a size-approximation algorithm which returns a u > a\ such that min{a : e > 8e(a)} > a u > ct m i n {e) 
and a m i n {e) > a\ > m&x{a : e < \e m i n (a)}. However, all these approximation algorithms still take high 
running time of 0(m 3 n 4 ) (the notation O hides poly log factors in m and n). 

Heuristics for tolerant-LCP. In practice, the tolerant-LCP is solved heuristically by using the geometric 
hashing and GHT algorithms for which rigorous analyses are only known for the exact-LCP. For example, 
the algorithms in [17, 31] are for tolerant-LCP but the analyses are for exact-LCP only. Because of its prac- 
tical performance, the exact version of GHT was carefully analyzed by Akutsu et al. [5], and a randomized 
version of the exact version of geometric hashing in 2D was given by Irani and Raghavan [26]. The tolerant 
version of GHT (and geometric hashing) is based on the corresponding exact version by replacing the exact 
matching with the approximate matching which requires a distance measure to compare the keys. We can 
no longer identify the optimal rigid motion by the maximum votes as in the exact case. Instead, the tolerant 
version of GHT clusters the rigid motions (which are points in a six-dimensional space) and heuristically 
approximates the optimal rigid motion by a rigid motion in the largest cluster. Thus besides not giving any 
guarantees about the solution, this heuristic requires clustering in six dimensions, which is computationally 
expensive. 

Other Related Work. There is some closely related work that aims at computing the minimum Hausdorff 
distance for PM (see, e.g., [13] and references therein). Also, the problems we are considering can be 
thought of as the point pattern matching problem under uniform distortion. Recently, there has been some 
work on point pattern matching under non-uniform distortion [28, 6]. 

Our results. There are three results in this paper. First, we introduce a new distance-approximation algo- 
rithm for tolerant-LCP algorithm, called Diheda (because our algorithm is based on Dihedral Angle 
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comparisons). 

Theorem 1.1 Let P, Q £ I 3 of size m and n, with m > n, and e > 0. Suppose that interpoint distances in 
P and in Q be > 2e (this is the condition for tolerant-LCP). DlHEDA (see Algorithm 1) finds a rigid motion 
fi and a subset I of Q such that 

• |/| > LCP(P,Q) and 

• dist(P, n{I)) < 4e 
in 0(m 3 n 3 log m) time. 

DlHEDA is simple and more efficient than the known distance-approximation algorithms (which are 
alignment-based) for tolerant-LCP. The running time of DlHEDA is 0(m 3 n 3 logm) in the worst case. For 
general input, we expect the algorithm to be much faster because it is simpler and more efficient than the 
previous heuristics that are known to be fast in practice. This is because our clustering step is simple (sorting 
linearly ordered data) while the clustering step in those heuristics requires clustering high-dimensional data. 

Second, based on a combinatorial observation, we improve the algorithms for exact-a-LCP by a linear 
factor for pose clustering or GHT and a quadratic factor for alignment or geometric hashing. This also 
corrects a mistake by Irani and Raghavan [26]. 

Finally, we achieve a similar speed-up for DlHEDA using a sampling approach based on expander graphs 
at the expense of approximation in the matched set size. We remark that this result is mainly of theoretical 
interest because of the large constant factor involved. Expander graphs have been used before in geometric 
optimization for fast deterministic algorithms [2, 27] ; however, the way we use these graphs appears to be 
new. Our results also hold when we extend the set of transformations to scaling; for simplicity we restrict 
ourselves to rigid motions in this paper. 

Outline. The paper is organized as follows. The rest of this section contains some preliminaries. In 
Section 2 we introduce our new distance-approximation algorithm for tolerant-LCP. In Section 3 we show 
how a simple deterministic sampling strategy based on the pigeonhole principle yields speed-ups for the 
exact-a-LCP algorithms. In Section 4 we show how to use expander graphs to further speed up the DlHEDA 
algorithm for tolerant-a-LCP at the expense of approximation in the matched set size. Section 5 is the 
conclusion. In the appendix, we recall and compare the existing four basic algorithms for exact-LCP: pose 
clustering, alignment, GHT and geometric hashing. 

Terminology and Notation. For a transformation p, denote by 1^ the set of points in p(Q) that are within 
distance e of some point in P. We call I M the matched set of p and say that p is an |7 M | -matching. We 
call the transformation p that maximizes the maximum matching transformation. A basis is a minimal 
(for containment relation) ordered tuple of points which is required to uniquely define a rigid motion. For 
example, in 2D every ordered pair is a basis; while in 3D, every non-collinear triplet is a basis. In Figure 1, 
a rigid motion in 3D is specified by mapping a basis (qi,q2, qs) to another basis (pi,P2,P3)- We call a key 
used to represent an ordered tuple S a rigid motion invariant key if it satisfies the following: (1) the key 
remains the same for all p(S) where p is any rigid motion, and (2) for any two ordered tuples S and S' with 
the same rigid motion invariant key there is a unique rigid motion p such that p(S) = S'. For example, 
as rigid motion preserves orientation and distances among points, given a non-degenerate triangle A, the 3 
side lengths of A together with the orientation (the sign of the determinant of the ordered triplet) form a 
rigid motion invariant key for A in M 3 . Henceforth, for simplicity of exposition, in the description of our 
algorithms we will omit the orientation part of the key. 
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Figure 1: In this example, the rigid motion is obtained by matching Cq = {q\, q%, . . . , #5} in Q to 
Cp = {pi,P2, ■ ■ ■ ,Ps} hi P- We have LCP(i- > , Q) = \Cp\ = \Cq\ = 5. The corresponding 5-matching 
transformation \l can be discovered by matching (^1,^2) to (pi,P2)> the rigid motions /ij that transform 
(qi, q2, q%) to (pi,P2,Pi) for i = 3, 4, 5 are all the same and thus \i = ^3 = ^4 = ^5 will get 3 votes, which 
is the maximum. 



2 DlHEDA 



In this section, we introduce a new distance-approximation algorithm, called DlHEDA , for tolerant-LCP. 
The algorithm is based on a simple geometric observation. It can be seen as an improvement of a known 
GHT-based heuristic such that the output has theoretical guarantees. 



2.1 Review of GHT 

First, we review the idea of the pair-based version of GHT for exact-LCP. See the appendix or [5, 31] for 
more details. For each congruent pair, say (pi,P2) m P an d (qi, 92) in Q, and for each of the remaining 
points p £ P and q 6 Q, if (qi,q2,q) is congruent to (pi,P2,p), compute the rigid motion fi that matches 
(<Zi> 92) q) to (pi,P2,p)- We then cast one vote for /x. The rigid motion that receives the maximum number 
of votes corresponds to the maximum matching transformation sought. See Figure 1 for an example. 



2.2 Comparable rigid motions by dihedral angles 

For the exact-LCP, one only needs to compare rigid motions by equality (for voting). For the tolerant-LCP, 
one needs to measure how close two rigid motions are. In IR 3 , each rigid motion can be described by 6 
parameters (3 for translations and 3 for rotations). How to define a distance measure between rigid motions? 
We will show below that the rigid motions considered in our algorithm are related to each other in a simple 
way that enables a natural notion of distance between the rigid motions. 



Observation. In the pair-based version of GHT as described above, the rigid motions to be compared have 
a special property: the rigid motions transform a common pair — they all match (qi, q^) to {p\,P2) in Fig- 
ure 1. Two such transformations no longer differ in all 6 parameters but differ in only one parameter. To see 
this, we first recall that a dihedral angle is the angle between two intersecting planes; see Figure 2 for an ex- 
ample. In general, we can decompose the rigid motion for matching (qi,q2, 93) to (pi,P2,P3) into two parts: 
first, we transform (gi, 92) to (pi,j)2) by a transformation fa; then we rotate the point ^1(^3) about P1P2 by 
an angle 9, where 9 is the dihedral angle between the planes (pi,P2,P3) and (fa(qi), fa(q2), fa{Q3))- This 
will bring 53 to coincide with p$. Thus, we have the following lemma: 
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. .. B - a ball around p with radius =£ 




\/ C - a circle 

Figure 2: The dihedral angle is the angle between planes formed by (pi,p2, q) and (pi,P2,p)- The rotation 
angles of transformations that rotate q about pip2 to within e of p form a subinterval of [0, 2n). 



Lemma 2.1 Let (pi,P2,P3) and (qi, q2, qz) be two congruent non-collinear triplets, and let <p\ be a rigid 
motion that takes qi to Pi for i = 1,2. Let <j)2 be the rotation about p\P2 by an angle 9, where 8 is the 
dihedral angle between the planes ipi,P2iPz) and (<Pi(qi), 4>i(Q2), 4>i(Q3))- Then the unique rigid motion 
that takes (pi,P2,P3) to (qi,q2, Q3) is equal to (j) 2 4>i- 

We now state another lemma that will be useful in the description and proof of correctness of DlHEDA . Let 
(pi ,p2,p) and q be four points as shown in Figure 2. Consider the rotations about P1P2 that take q to within 
e of p. The rotation angles of these transformations form a subinterval of [0, 2ir). This is because a circle 
C (corresponding to the trajectory of p) intersects with the sphere B (around p with radius e) at at most 
two points (corresponding to a subinterval of [0, 2n)), as shown in Figure 2. That is, we have the following 
lemma: 

Lemma 2.2 Let pi,p2,p, q £M. 3 be four points (not necessarily non-collinear), then the rotation angles of 
transformations that rotate q about p\P2 to within e of p form a subinterval of [0, 2tt). 

2.3 Approximating the optimal rigid motion by the "diametric" rigid motion 

For a point set S C R 3 , we call a pair of points {p,q} € S 2 diameter-pair if \\p — q\\ = diameter(S'). Arigid 
motion of Q that takes q\ to p\ and 52 on the line p\P2 and closest possible to P2 is called a (pi,P2, Qi, Q2)- 
rigid motion. Based on an idea similar to the one behind Lemma 2.4 in Goodrich et al. [19], we have the 
following lemma: 

Lemma 2.3 Let p be a rigid motion such that each point of p(S), where S C Q, is within distance e of a 
point in P. Let 92} be a diameter-pair of S. Let pi G P be the closest point to fi(qi)for i = 1,2. Then 
we have a (pi,P2, Qi,Q2)-rigid motion p! of Q such that each point of p,'(S) is within 4e of a point in P. 

Proof Sketch. Translate p(q\) to p\, this translation shifts each point by at most e. Next, rotate about p\ such 
that p(q2) is closest to P2 (which implies //(<?i), At'fe) and P2 are collinear). Since {^i, 52} is a diameter- 
pair, this rotation moves each point by at most 2e. Thus, each point is at most e + e + 2e = 4e from its 
matched point. ■ 

2.4 Approximation algorithm for tolerant-LCP 

We first describe the idea of our algorithm DlHEDA . Input is two point sets in M 3 , P = {pi, . . . ,p m } and 
Q = {^i i • • • > Qn} with m > n, and e > 0. Suppose that the optimal rigid motion was achieved by 
matching a set 7 W) = {gi, q 2 , . . . , q k } Q Q to J M = {pi,p 2 , ■ ■ ■ ,Pk} Q P- WLOG, assume that {q 1: q 2 } 
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is the diameter pair of Then by Lemma 2.3, there exists a (pi,P2, Qi, Q2 ) -rigid motion /i of Q such that 
fj,(I^ ) is within 4e of a point in P. Since we do not know the matched set, we do not know a diameter-pair 
for the matched set either. Therefore, we exhaustively go through each possible pair. Namely, for each pair 
((7i> 92) £ Q and each pair (pi,P2) £ P, if they are approximately congruent then we find a (pi,P2, 9i> (Z2)- 
rigid motion /x of Q that matches as many remaining points as possible. Note that (pi,p2,(?i,(Z2)-rigid 
motions are determined up to a rotation about the line p±p2- By Lemma 2.2, the rotation angles that bring 
fi(qi) to within 4e of pi form a subinterval of [0, 2ir). And the number of non-empty intersection subintervals 
corresponds to the size of the matched set. Thus, to find fi, for each pair (p,q) € P\ {pi , P2 } x Q \ {qi , 12 }, 
we compute the dihedral angle interval according to Lemma 2.2. The rigid motion /j sought corresponds to 
an angle <p that lies in the maximum number of dihedral intervals. The details of the algorithm are described 
in Algorithm 1 . 

Time Complexity. For each triplet in Q, using kd-tree for range query, it takes 0(m 3 '( 1_ 3) + m 3 + 
m 3 log m 2 ) = 0(m 3 log to) for lines 1 1-20. For each pair (q\ , 52) and (p\ , P2), we spend time 0(mn) to 
find the subintervals for the dihedral angles, and time 0{mn log m) to sort these subintervals and do the scan 
to find an angle that lies in the maximum number of subintervals. Thus the total time is 0(m 3 n 3 log m). 

3 Improvement by pigeonhole principle 

In this section we show how a simple deterministic sampling strategy based on the pigeonhole principle 
yields speed-ups for the four basic algorithms for exact-a-LCP. Specifically, we get a linear speed-up for 
pose clustering and GHT, and quadratic speed-up for alignment and geometric hashing. It appears to have 
been erroneously concluded previously that no such improvements were possible deterministically [26]. 

In pose clustering or GHT, suppose we know a pair (qi , (72 ) in Q that is in the sought matched set, then the 
transformation sought will be the one receiving the maximum number of votes among the transformations 
computed for (q\, 52)- Thus if we have chosen apair (qi, 92) that lies in the matched set, then the maximum 
matching transformation will be found. We are interested in the question "can we find a pair in the matched 
set without exhaustive enumeration"? The answer is yes: we only need to try a linear number of pairs 
(<Zi> 92) to find the maximum matching transformation or conclude that there is none that matches at least g 
points. 

We are given a set Q = {qi, . . . , q n }, and let I C Q be an unknown set of size > g for some constant 
a > 1. We need to discover a pair (p, q) with p,q £ I by using queries of the following type. A query 
consist of a pair (a, b) with a, b G Q. If we have a,b G /, the answer to the query is YES, otherwise the 
answer is NO. Thus our goal is to devise a deterministic query scheme such that as few queries are needed 
as possible in the worst case (over the choice of /) before a query is answered YES. Similarly, one can ask 
the question about querying triplets to discover a triplet entirely in /. 

Theorem 3.1 For an unknown set I C Q with \I\ > ^ and \Q\ = n using queries as described above, 

(1) it suffices to query 0(an) pairs to discover a pair in I; 

(2) it suffices to query 0(a 2 n) triplets to discover a triplet in I. 

Proof. The proof is based on the pigeonhole principle. To prove (1), we assume for simplicity that a and ^ 
are both integers. Partition the set Q into ^ subsets of size a each. Since the size of / is more than ^, by the 
pigeonhole principle, there is a pair of points in / that lies in one of the above chosen subsets. Thus querying 
all pairs in these subsets will discover /. This gives that - (2) ~ an queries are sufficient to discover /. 

Similarly, to prove (2), partition Q into ^- subsets Pi, ... , Pjl. of size 2a each (we assume, as before, 
that 2a and ^ are both integers). Now we test all triplets that lie in the Pi's. Any set I C Q that intersects 
with each of the Pj's in at most 2 points has size < ^. Hence if |/| > ^ then it must intersect with one 
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Algorithm 1 DlHEDA 



1: procedure PREPROCESSING 

2: for each pair (p\ , p 2 ) of P do 

3: Compute and insert the key of | \p1P2 1 1 into a dictionary V\ ; 

4: end for 

5: for each triplet (pi , p 2 , P3) of P do 

6: Compute and insert the n'gj'<i motion invariant key for (pi , p 2 , P3) into a dictionary D 2 ; 

1: end for 

8: end procedure 

9: procedure Recognition 

10: for each pair (q\ , g 2 ) of (^) do o This can be reduced by the edge set of an expander of Q. 

11: if [ 1 1 <?i <?2 1 1 — 2e, ||<zi<72|| + 2e] exists in V\ then 

12: Initialize an empty dictionary P3 of pairs; 

13: for each remaining point q G Q do 

14: Compute and search the range [ 1 1 <?i <?2 1 1 — 2e, H^i^H + 2e] x [||<7<7i|| — 2e, | l^gi 1 1 + 2e] x 
OI992II - 2e, ||q^ 2 || + 2e] of {qi,q2,q) in £> 2 ; > e.g. using a kd-tree. 

15: for each entry (pi ,p 2 ,p) found do 

16: If {pi,P2) exists in P3, increase its vote; otherwise insert (pi,p 2 ) into P3 with one 

vote; 

17: Append the matched pair (q,p) to the list associated with (pi,p 2 )', 

18: end for 

19: end for 

20: end if > Compute the maximum transformation that matches (qi , g 2 ) to (p\ ,p 2 ). 

21: for each pair (pi , p 2 ) in the dictionary D3 do 

22: Compute a transformation (f> that brings gi to p\ and g 2 closest to p 2 ; 

23: For each matched pair (q, p) of the associated list of (pi,p 2 ), compute an interval of dihedral 

angles such that <j)(q) is within 4e of p; 

24: Sort all the intervals of dihedral angles; and find a dihedral angle tp that occurs in the largest 

number V of intervals; 

25: Compute the transformation /x by the composition of <fi and the rotation about PYP2 by angle 
ip; > n brings V + 2 points of Q to within 4e of some matched points in P. 

26: Keep the maximum matched set size and the corresponding transformation; 

27: end for 

28: end for 

29: end procedure 
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of the sets above in at least 3 points. Thus testing the triplets from the Pi's is sufficient to discover /. The 
number of triplets tested is ^ ( 2 ^) ~ a 2 n. ■ 

Remark: It can be shown that the schemes in the proof above are the best possible in requiring the 
smallest number of queries (up to constant factors). 

In alignment and geometric hashing algorithms if we have chosen a triplet (91 , 92 , 93) from the maximum 
matching set / C Q then we will discover /. The question, as before, is how many triplets in Q need to be 
queried to discover a set I of size > ^. By Theorem 3.1 (2), we only need to query 0(a 2 n) triplets. Thus 
the running times of both alignment and geometric hashing are improved by a factor of 6(n 2 ). 

See Table 1 for the time complexity comparison of deterministic algorithms for exact-a-LCP in R 3 . 

Finally, our approximation algorithm for tolerant-LCP adapts naturally for exact-a-LCP with pigeonhole 
sampling. We analyze the running time of our algorithm for exact-a-LCP with the pigeonhole sampling of 
pairs. In the exact case, each exact matched pair of points (q,p) corresponds to a single dihedral angle. We 
thus find the dihedral angle that occurs the maximum number of times by sorting all the dihedral angles. For 
a fixed pair (91, 92) and a point q in Q the number of triplets ((pi,P2),P3) in P that match ((91, 92), 93) is 
bounded above by 3H 2 (m), where ^(m) is the maximum possible number of the congruent triangles in a 
point set of size m in R 3 . Total time spent for pair (91,92) then is 0(nPt2). Since we use 0(an) pairs, the 

overall running time is 0{an 2 H2). Agarwal and Sharir [1] show that H 2 (m) < mzg(m), where g(m) is a 
very slowly growing function of m of inverse-Ackermann type. 



Algorithm 


Original running time 


Improved running time 


Pose Clustering (e.g. [31]) 


0(rn i n i S(m)) 


0(m i n 2 S(m)) 


Alignment (e.g. [5]) 


0(m A + m\ A,2 (m, n))S(m) 


0(m A n 2 S(m)) 


GHT (e.g. [5]) 


0{m 6 S{m) + A 3 > 2 (m, n)S(X s ' 2 (m, n))) 


0(m A S(m) + n 2 H 2 (m)) 


Geometric hashing (e.g. [30]) 


0{m A s\m) +n i H 3 (m)) 


0(m 4 S(m) +n 2 H 3 (m)) 


This paper 




o{m 6 s{m)+ri 2 H 2 {m)) 



Table 1: Time complexity comparison of deterministic algorithms for exact-a-LCP in M 3 . S(x) is the query 
time for the dictionary of size x, which can be taken to be O(logx) or smaller; H2(m) is the maximum 
number of triangles spanned by m points in M 3 that are congruent to a given triangle, we have H2(m) < 
m3 g(m), where g(m) is a very slowly growing inverse-Ackermann type function of m [1], and can be 
treated as constant for all practical purposes; H^{m) is the maximum number of tetrahedrons spanned by 
n points in M 3 that are congruent to a given tetrahedron, we have H^(m) = 0(m 2+s ) for any 5 > [1]; 
A 3 ' 2 (m,ra) = 0(min{m L8 n 3 ,m L95 n 2 - 68 + m L8875 ra 2 - 8 })/5/. 

As is often the case for algorithms for LCP, analysis involves determining quantities such as H2(m), 
which is a difficult problem. In the above table we have tried to give references for the first four algorithms 
including the tightest analyses rather than the original sources. Note that our algorithm is simpler than the 
others in the first column which involve checking for congruent simplices in a dictionary. 

4 Expander-based sampling 

While for the exact-a-LCP the simple pigeonhole sampling served us well, for the tolerant-a-LCP we do 
not know any such simple scheme for choosing pairs. The reason is that now we not only need to guarantee 
that each large set contain some sampled pairs, but also that each large set contain a sampled pair with large 
length (diameter-pair) as needed for the application of Lemma 2.3 in the DlHEDA algorithm. Our approach 
is based on expander graphs (see, e.g., [3]). Informally, expander graphs have linear number of edges but the 
edges are "well-spread" in the sense that there is an edge between any two sufficiently large disjoint subsets 
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of vertices. Let G be an expander graph with Q as its vertex set. We show that for each S C Q, if \S\ is not 
too small, then there is an edge (u, v) in G such that [u, v) G S 2 and | \uv\ | approximates the diameter of S. 

By choosing the pairs for the DlHEDA algorithm from the edge set of G (the rest of the algorithm is same 
as before), we obtain a bicriteria - distance and size - approximation algorithm as stated in Theorem 4.4 
below. We first give a few definitions and recall a result about expander graphs that we will need to prove 
the correctness of our algorithm. 

Definition 4.1 Let S be a finite set of points of W for r > 1, and let < k < n. Define diameter^, k) = 
min T .| T | =fc diameter^ \ T). 

That is, diameter^, k) is the minimum of the diameter of the sets obtained by deleting k points from S. 
Clearly, diameter^, 0) = diameter(S'). 

Let U and V be two disjoint subsets of vertices of a graph G. Denote by e(U, V) the set of edges in G 
with one end in U and the other in V. We will make use of the following well-known theorem about the 
eigenvalues of graphs (see, e.g. [29], for the proof and related background). 

Theorem 4.2 Let G be a d-regular graph on n vertices. Let d = \\ > A2 > • • • > A n be the eigenvalues of 
the adjacency matrix of G. Denote A = max2<j< n |Aj|. Then for every two disjoint subsets U,W C V, 

n 

Corollary 4.3 Let U, W C V be two disjoint sets with \U\ = \W\ > Then G has an edge in U xW. 

Proof. It follows from (1) that if > A^/jfTp 7 ! then \e(U,W)\ > 0, and since \e(U,W)\ is 

integral, \e(U, W)\ > 1. But the above condition is clearly true if we take U and W as in the statement of 
Corollary 4.3. ■ 

There are efficient constructions of graph families known with A < 2\fd (see, e.g., [3]). Let us call such 
graphs good expander graphs. We can now state our main result for this section. 

Theorem 4.4 For an a-LCP instance (P, Q) with LCP(P, Q) > ^, the DlHEDA algorithm with expander- 
based sampling using a good expander graph of degree d > 2500a 2 finds a rigid motion /i in time 
0(m 3 n 2 logm) such that there is a subset I satisfying the following criteria: 

(1) size-approximation criterion: \I\ > LCP(P,Q) — 

(2) distance-approximation criterion: each point of p,(I) is within distance Qefrom a point in P. 

Thus by choosing d large enough we can get as good size-approximation as desired. The constants in the 
above theorem have been chosen for simplicity of the proof and can be improved slightly. 

For the proof we first need a lemma showing that choosing the query pairs from a graph with small X(G) 
(the second largest eigenvalue of G) gives a long (in a well-defined sense) edge in every not too small subset 
of vertices. 

Lemma 4.5 Let G be a d-regular graph with vertex set Q C M 3 , and \Q\ = n. Let S C Q be such that 
\S\ > 25A ^ Tt . Then there is an edge {s±, S2} G E(G) n S 2 such that ||siS2|| > diameter ^' — 4 — ™)_ 

Proof. For a positive constant c to be chosen later, remove cn pairs from S as follows. First remove a 
diameter pair, then from the remaining points remove a diameter pair, and so on. Let T be the set of points in 
the removed pairs and T p the set of removed pairs. The remaining set S\T has diameter > diameter(S, 2cn) 
by the definition of diameter^, 2cn), and hence each of the removed pairs has length > diameter^, 2cn). 
For B,C C 5 let ||B,C|| = min 6eB)CeC ||6c||. 



< \V\U\\W\. (1) 
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Claim 1 The set T defined above can be partitioned into three sets B, C, E, such that \B\, \C\ > ^P, and 
\\B C\\ > d ' ameter > 2cn ) 

Proof. Fix a Cartesian coordinate system and consider the projections of the pairs in T p on the x-, y- 
,and z-axes. It is easy to see that for at least one of these axes, at least ^ pairs have projections of length 
> dBmet ^ S|2m ) , Suppose without loss of generality that this is the case for the x-axis, and denote the set of 
projections of pairs on the x-axis with length > diame W 5 ' 2cw ) by T%, and the set of points in the pairs in T% 

by T x . We have \T X \ > 2cn/3. Now consider a sliding window W on the x-axis of length d ' amete ^( s > 2cn ) ; 
initially at — oo, and slide it to +00. At any position of W, each pair in T% has at most 1 point in W, as the 
length of any pair is more than the length of W. Thus at any position, W contains < \T X \ = \T x \/2 points. 
It is now easy to see by a standard continuity argument that there is a position of W, call it W, where there 
are > > m. points of T x both to the left and to the right of W. 

Now, B is defined to be the set of points in T whose projection is in T x and is to the left of W; similarly 
C is the set of points in T whose projection is in T x and is to the right of W. Clearly any two points, one 
from B and the other from C, are d ' amete ^ S ' 2cn ) -apart. ■ 

Coming back to the proof of Lemma 4.5, the property that we need from the query-graph is that for any two 
disjoint sets B, C C S of size 5\S\, where S is a small positive constant, the query-graph should have an 
edge inBxC. 

By Corollary 4.3 if \B\ > f > ?f, and \C\ > f > ?f, that is, ifof, then G has an edge inBxC. 
Taking c = 12 f X completes the proof of Lemma 4.5. ■ 

Proof of Theorem 4.4. If we take G to be a good expander graph then Lemma 4.5 gives that G has an 

diameter(S,^=n) 

edge of length > ^ — ^ — L e t g a lso be a solution to tolerant-LCP for input (P,Q) with error 

diameter^, ^=n) 

parameter e > 0. We have that one of the sampled pairs has length at least ^ — — ■ Thus applying 

an appropriate variant (replacing the diameter pair by the sampled pair with large length as guaranteed by 
Lemma 4.5) of Lemma 2.3, we get a rigid motion /i such that there is a subset / satisfying the following: 

(1) \I\ > \S\ - for any d > 2500a 2 ; 

(2) Each point of / is within 6e(= e + e + 4e) of a point in M. ■ 



5 Discussion 

We have presented a new practical algorithm for point pattern matching. Our DlHEDA algorithm is the 
fastest known distance-approximation algorithm for tolerant-LCP, and is simple compared to other known 
distance-approximation algorithms and heuristics which involve 6-dimensional clustering. Our analysis of 
DlHEDA is not tight, and perhaps better bounds can be obtained if the interpoint distance is greater than e 
by a sufficiently large constant factor. 

Our technique of pigeonhole sampling yields speed-ups for all four popular algorithms and also the 
fastest known deterministic algorithm for the exact-LCP. Again, our algorithms are simpler than the previous 
best algorithms. Akutsu et al. [5] give a tighter analysis for GHT in terms of the function A 3 ' 2 (m, n). Our 
analysis of DlHEDA (and GHT) with pigeonhole sampling was based on H2(m). Presumably, a better 
analysis similar to the idea in [5] is possible. 

Point pattern matching is of fundamental importance for computer vision and structural bioinformatics. 
Indeed, this investigation stemmed from research in structural bioinformatics. Current software, which uses 
either geometric hashing or generalized Hough transform, can immediately benefit from this work. We have 



11 



implemented a randomized version of DlHEDA for molecular common substructure detection and the results 
were reported in [16]. 

Acknowledgment. We thank S. Muthukrishnan and Ali Shokoufandeh for the helpful comments and 
advice. 
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Appendix 



A Voting Algorithms for Exact-LCP 

In this appendix, we review and compare four popular algorithms for exact-LCP: pose clustering, alignment, 
generalized Hough transform(GWY), and geometric hashing. These algorithms are all based on a voting idea 
and are sometimes confused in the literature. Please see Algorithms 2, 3, 4, 5) for a full description of the 
algorithms in their generic form independent of the search data structure used. In particular, geometric hash- 
ing algorithms need not use a hash-table as a search data structure. We describe all the algorithms in terms 
of a dictionary of objects (which are either transformations or a set of points and can be ordered lexico- 
graphically). Denote the query time for this dictionary by S(x) + O(k) where x is the size of the dictionary, 
and k is the size of the output depending on the query. For example, if the dictionary is implemented by a 
search tree we have S(x) = O(logx). 

Pose clustering and alignment are the basic methods. GHT and geometric hashing can be regarded as 
their respective efficient implementations. Efficiency is achieved by preprocessing of the point sets using 
their rigid motion invariant keys which speeds-up the searches. 

In pose clustering, for each pair of triplets (qi,q2, qs) € Q and (pi,P2,Ps) £ P, we check if they are 
congruent. If they are then we compute the rigid motion jj, such that fi(qi,q2, q%) = (pi,P2,Ps)- We then 
cast one vote for /j,. The rigid motion which receives the maximum number of votes corresponds to the 
maximum matching transformation sought. The running time of pose clustering is 0(m 3 n 3 S(m 3 n 3 )) as 
the size of the dictionary of transformations can be as large as 0(m 3 n 3 ). 

In alignment, for each pair of triplets (</i,<?2, #3) £ Q and (pi,P2,P3) £ P we check if they are congru- 
ent. If they are then we compute the rigid motion /j, such that n(q±, q2, 93) = {pi,P2,Ps)- Then we count the 
number of points in fi(Q) that coincide with points in P. This number gives the number of votes the rigid 
motion fi gets. The rigid motion which receives the maximum number of votes corresponds to the maximum 
matching transformation sought. The running time is 0(m 3 n 4 S(m)). 

The difference between pose clustering and alignment is the voting space: in pose clustering voting is 
done for transformations while in alignment it is for bases (triplets of points). In both pose clustering and 
alignment algorithms, each possible triplet in Q is compared with each possible triplet in P. However, by 
representing each triplet with its rigid motion invariant key, only triplets with the same key (rigid motion 
invariant) are needed to be compared. This provides an efficient implementation. For example, the GHT 
algorithm is an efficient implementation of pose clustering. Here we preprocess P by storing the triplets 
of points with the rigid motion invariant keys in a dictionary. Now for each triplet (q±, q2, qs) in Q we 
find congruent triplets in P by searching for the rigid motion invariant key for (q±, q2, #3). The rest of 
the algorithm is the same as pose clustering. Similarly the geometric hashing algorithm is an efficient 
implementation of the alignment method. 

GHT is faster than geometric hashing, however geometric hashing has the advantage that algorithm can 
stop as soon as it has found a good match. Depending on the application this gives geometric hashing 
advantage over GHT. 

As observed by Olson [31] and Akutsu et al. [5], pose clustering and GHT can be further improved. This 
is because a /c-matching transformation can be identified by matching (k — 2) bases which match a common 
pair. We call this version of the generalized Hough transform the pair-based version; it is described below 
in Algorithm 6. Although the worst case time complexity of the pair-based version and the original version 
are the same, this will serve as a basis for our new scheme, called DlHEDA . The pair-based version also 
allows efficient random sampling of pairs [31,5]. 
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Algorithm 2 Pose Clustering 



procedure Pose Clustering^, Q) 

Initialize an empty dictionary V of rigid motions; 
for each triplet (q±, q 2 , 93) of Q do 
for each triplet (pi,P2,P3) of P do 

if (qi, q 2 , qs) is congruent to (pi,P2,Ps), then 

Compute the rigid motion fi which matches (qi,q2, qs) to (pi,P2,P3) ; 
Search /i in the dictionary V; 

If found, increase the votes of fi; otherwise insert fi with one vote, 
end if 
end for 
end for 

Return the maximum vote rigid motion in V; 
end procedure 



Algorithm 3 Alignment 
1: procedure Alignment(P, Q) 



2: for each triplet (qi,q 2 , #3) of Q do 
3: for each triplet (pi , P2, P3) of P do 

4: If (qi, q2, (fe) is congruent to (pi,P2,P3), compute the rigid motion \i\ 

5: Vote = 0; > Vote is a local counter for the transformation fi. 

6: for each remaining point q G Q and p G P do 

7: If /[/(g) = p, then increase Vote by 1; 

8: end for 

9: Keep the maximum vote and its associated transformation; 

10: end for 

ll: end for 

12: Return the maximum vote transformation. 



13: end procedure 
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Algorithm 4 The original version of generalized Hough transform. 



1: procedure PREPROCESSING 

2: for each triplet (jp\ , p2 , P3) of P do 

3: Compute and insert the rigid motion invariant key for (pi ,P2,Ps) into a dictionary T>\ ; 

4: end for 

5: end procedure 

6: procedure RECOGNITION 

7: Initialize an empty dictionary T>2 of rigid motions; 
8: for each triplet (^i, ^2, 93) of Q do 

9: Compute and search the rigid motion invariant key for (qi,q2, 93) in the dictionary V\\ 

10: for each entry (pi,P2,P3) found, do 

11: Compute the rigid motion /j, which matches (q\,q2, 93) to (pi,P2,P3)', 

12: Search /x in the dictionary D2; 

13: If found, increase the votes of /x; otherwise insert /x with one vote into P2; 

14: end for 

15: end for 

16: Return the maximum vote rigid motion in T>2\ 

17: end procedure 



Algorithm 5 Geometric Hashing 

1: procedure PREPROCESSING 

2: for each triplet (pi , P2 , ps) of P do 

3: for each of the remaining point p of P do 

4: Compute and insert the rigid motion invariant key for {(pi, P2, P3), p} into a dictionary V\ ; 

5: end for 

6: end for 

7: end procedure 

8: procedure Recognition 

9: for each triplet (gi, q2, 53) of Q do 
10: Build an empty dictionary T>2 (of triplets of P); 

11: for each of the remaining point q of Q do 

12: Compute and search the rigid motion invariant key for {(gi, q2, 53), g} in the dictionary Pi; 

1 3 : for each entry {(pi,P2,P3), p} found do 

14: If (pi,P2,P3) exists in T>2, then increase its vote by one; otherwise insert (^1,^2,^3) into 

T>2 with vote one. 
15: end for 

16: end for 

17: Keep the maximum vote and compute the corresponding transformation from its associated 

triplet; 
18: end for 

19: end procedure 
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Algorithm 6 The pair-based version of generalized Hough transform. 



1: procedure PREPROCESSING 

2: for each pair (p\,p2) of P do 

3: for each remaining point p of P do 

4: Compute and insert the rigid motion invariant key for {{pi,P2),p} into a dictionary V; 

5: end for 

6: end for 

7: end procedure 

8: procedure Recognition 

9: for each pair (q\ , (ft) of Q do 

10: Initialize an empty dictionary T>i of rigid motions; 

11: for each remaining point q of Q do 

12: Compute and search the rigid motion invariant key for {(qi, (ft), 9} in the dictionary 79; 

13: for each entry {{pi,P2),p} found, do 

14: Compute the rigid motion p which matches {(qi, (ft), 9} to {(pi,P2),p}- 

15: Search /x in the dictionary T>2\ 

16: If found, increase the votes of p\ otherwise insert p with one vote into T>2\ 

17: end for 

18: end for 

19: Keep the rigid motion for (gi, 92) that receives the maximum number of votes. 

20: end for 

21: Return the rigid motion that receives the maximum number of votes among all pairs. 

22: end procedure 
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